
#This module contains user defined functions (casadi) that which relate volumetric moisture, capillary capacity, and the hydraulic conductivity with pressure head (h)

from casadi import *

# thetaR = 0.078
# thetaS = 0.43
# alpha = 3.6
# n = 1.56
# m = 1-1/n
# Ks = 0.00000288889
# neta = 0.5
# Ss = 0.00001

thetaR = 0.065  # 0.078
thetaS = 0.41  # 0.43
alpha = 7.5  # 3.6
n = 1.89  # 1.56
m = 1-1/n
Ks = 4.4208e-2/3600  # 0.00000288889
neta  = 0.5
Ss = 0.00001


def volumetricMoisture(head):
    Se = (1 + fabs(head*alpha)**n)**(-m)
    volMoisture = thetaR + (thetaS - thetaR)*Se 
    return volMoisture
	#Dependence of volumetric moisture content on the pressure head
    
# 	 Se = if_else(head >= 0.,1.,(1 + fabs(head*alpha)**n)**(-m))


def capillaryCapacity(head):
	#Dependence of the capillary capacity on the pressure head
    Se = (1 + fabs(head*alpha)**n)**(-m)
    dSedh = alpha*m/(1 - m)*Se**(1/m)*(1-Se**(1/m))**m
    capCapacity = Se*Ss+(thetaS-thetaR)*dSedh
    return capCapacity
	


def hydraulicConductivity(head):
	#Dependence of unsaturated hydraulic conductivity on the pressure head
    Se = (1 + fabs(head*alpha)**n)**(-m)
    hyrdConductivity = Ks*Se**neta*(1-(1-Se**(1/m))**m)**2
    return hyrdConductivity
	
	
# 	Se=if_else(head>=0.,1.,(1+fabs(head*alpha)**n)**(-m))
def pressureHead(volMoisture):
	#Dependence of pressure head on the volumetric moisture content
	presHead = (((((volMoisture - thetaR) / (thetaS - thetaR + 1.e-20) + 1.e-20) ** (1. / (-(1-1/(n+1.e-20)) + 1.e-20))
              - 1) + 1.e-20) ** (1. / (n + 1.e-20))) / (-alpha + 1.e-20)
	return presHead


def hydraulicConductivityApprox(head_1, head_2):
	#The approximation of the hydraulic conductivity at the center of the compartment
	hydCondApprox = 0.5*(hydraulicConductivity(head_1) + hydraulicConductivity(head_2))
	return hydCondApprox
